### When Democracy Brings Insecurity
### Hans Lueders
### Eurobarometer data
### 18 April 2023




# Preparation -------------------------------------------------------------

### packages
library(ggplot2)
library(stargazer)
library(lfe)
library(xtable)
library(readxl)
library(tm)
library(stringr)
library(SnowballC)
library(rgdal)
library(MASS)



### set working directory
path <- "" # define path here
setwd(path)

### empty environment
rm(list=ls())


### load data
load("Data/WhenDemocracyBringsInsecurity_HansLueders_Eurobarometer.RData")





# Table 9 -----------------------------------------------------------------

### Regressions
# satisfaction with country's current economic situation
m1 <- felm(as.numeric(demsatis) ~ as.numeric(econsit_country)*as.factor(region) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)
m2 <- felm(as.numeric(demsatis) ~ as.numeric(econsit_country)*as.factor(region)*born_before_1975 + 
             as.numeric(econsit_country)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# satisfaction with own current job situation
m3 <- felm(as.numeric(demsatis) ~ as.numeric(jobsit_personal)*as.factor(region) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)
m4 <- felm(as.numeric(demsatis) ~ as.numeric(jobsit_personal)*as.factor(region)*born_before_1975 + 
             as.numeric(jobsit_personal)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# expectations about country's future economic situation
m5 <- felm(as.numeric(demsatis) ~ as.numeric(expect_econsit)*as.factor(region) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)
m6 <- felm(as.numeric(demsatis) ~ as.numeric(expect_econsit)*as.factor(region)*born_before_1975 + 
             as.numeric(expect_econsit)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# expectations about own future job situation
m7 <- felm(as.numeric(demsatis) ~ as.numeric(expect_jobsit)*as.factor(region) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)
m8 <- felm(as.numeric(demsatis) ~ as.numeric(expect_jobsit)*as.factor(region)*born_before_1975 + 
             as.numeric(expect_jobsit)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)


### Stargazer output
stargazer(m1,m2,m3,m4,m5,m6,m7,m8,
          omit = c("female", "birthcohort", "agegroup", "employment", "marstat", "age_education"),
          omit.stat = c("rsq", "f", "ser"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(8),
          column.labels = c("(1=not at all satisfied; 4=very satisfied)"),          
          covariate.labels=c("Country's current economic situation",
                             "Own current job situation",
                             "Country's future economic situation",
                             "Own future job situation",
                             "Eastern Europe",
                             "Born before 1975",
                             "Country's current economic situation x Eastern Europe",
                             "Country's current economic situation x Born before 1975",
                             "Eastern Europe x Born before 1975",
                             "Country's current economic situation x Eastern Europe x Born before 1975",
                             "Own current job situation x Eastern Europe x Born before 1975",
                             "Own current job situation x Eastern Europe",
                             "Own current job situation x Born before 1975",
                             "Country's future economic situation x Eastern Europe x Born before 1975",
                             "Country's future economic situation x Eastern Europe",
                             "Country's future economic situation x Born before 1975",
                             "Own future job situation x Eastern Europe x Born before 1975",
                             "Own future job situation x Eastern Europe",
                             "Own future job situation x Born before 1975"), 
          digits = 3,
          add.lines = list(
            c("Controls included?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Country-FE?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Year-FE?", "YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")))





# Table 10 ----------------------------------------------------------------

### create new region variable
d$region_UER <- NA
d$region_UER[d$region == "Western/Northern Europe"] <- 0
d$region_UER[d$country %in% c("CZ - Czech Republic", "EE - Estonia", "HU - Hungary", "RO - Romania", "SI - Slovenia")] <- 1
d$region_UER[d$country %in% c("BG - Bulgaria", "LV - Latvia", "LT - Lithuania", "PL - Poland", "SK - Slovakia")] <- 2
d$region_UER <- factor(d$region_UER, levels=0:2,
                       labels = c("Western/Northern Europe", "Low UER", "High UER"))

### Regressions
# satisfaction with country's current economic situation
m1 <- felm(as.numeric(demsatis) ~ as.numeric(econsit_country)*as.factor(region_UER) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d, 
           subset = region=="Eastern Europe")

# satisfaction with own current job situation
m2 <- felm(as.numeric(demsatis) ~ as.numeric(jobsit_personal)*as.factor(region_UER) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d, 
           subset = region=="Eastern Europe")

# expectations about country's future economic situation
m3 <- felm(as.numeric(demsatis) ~ as.numeric(expect_econsit)*as.factor(region_UER) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d, 
           subset = region=="Eastern Europe")

# expectations about own future job situation
m4 <- felm(as.numeric(demsatis) ~ as.numeric(expect_jobsit)*as.factor(region_UER) + 
             as.factor(birthcohort) + as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d, 
           subset = region=="Eastern Europe")


### Stargazer output
stargazer(m1,m2,m3,m4,
          omit = c("female", "birthcohort", "agegroup", "employment", "marstat", "age_education"),
          omit.stat = c("rsq", "f", "ser"),
          star.cutoffs = c(0.05,0.01,0.001),          
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(8),
          column.labels = c("(1=not at all satisfied; 4=very satisfied)"),          
          covariate.labels=c("Country's current economic situation",
                             "Own current job situation",
                             "Country's future economic situation",
                             "Own future job situation",
                             "High UER",
                             "Country's current economic situation x High UER",
                             "Own current job situation x High UER",
                             "Country's future economic situation x High UER",
                             "Own future job situation x High UER"), 
          digits = 3,
          add.lines = list(
            c("Controls included?", "YES", "YES", "YES", "YES"),
            c("Country-FE?", "YES", "YES", "YES", "YES"),
            c("Year-FE?", "YES", "YES", "YES", "YES")))





# Table A24 ---------------------------------------------------------------

### write function
desc <- function(dv){
  # define DV
  d$dv <- as.numeric(d[,dv])
  
  # get mean
  mean <- mean(d$dv, na.rm=T)
  
  # get SD
  sd <- sd(d$dv, na.rm=T)
  
  # output
  return(cbind(mean,sd))
}

### create dummies for region and birth cohort
# region
d$region_CEE <- NA
d$region_CEE[d$region == "Western/Northern Europe"] <- 0
d$region_CEE[d$region == "Eastern Europe"] <- 1

d$region_West <- NA
d$region_West[d$region == "Western/Northern Europe"] <- 1
d$region_West[d$region == "Eastern Europe"] <- 0

# birth cohort
d$birth5year_1990 <- NA
d$birth5year_1990[d$birth5year != "1990+" & !is.na(d$birth5year)] <- 0
d$birth5year_1990[d$birth5year == "1990+" & !is.na(d$birth5year)] <- 1

d$birth5year_8589 <- NA
d$birth5year_8589[d$birth5year != "1985-89" & !is.na(d$birth5year)] <- 0
d$birth5year_8589[d$birth5year == "1985-89" & !is.na(d$birth5year)] <- 1

d$birth5year_8084 <- NA
d$birth5year_8084[d$birth5year != "1980-84" & !is.na(d$birth5year)] <- 0
d$birth5year_8084[d$birth5year == "1980-84" & !is.na(d$birth5year)] <- 1

d$birth5year_7579 <- NA
d$birth5year_7579[d$birth5year != "1975-79" & !is.na(d$birth5year)] <- 0
d$birth5year_7579[d$birth5year == "1975-79" & !is.na(d$birth5year)] <- 1

d$birth5year_7074 <- NA
d$birth5year_7074[d$birth5year != "1970-74" & !is.na(d$birth5year)] <- 0
d$birth5year_7074[d$birth5year == "1970-74" & !is.na(d$birth5year)] <- 1

d$birth5year_6569 <- NA
d$birth5year_6569[d$birth5year != "1965-69" & !is.na(d$birth5year)] <- 0
d$birth5year_6569[d$birth5year == "1965-69" & !is.na(d$birth5year)] <- 1

d$birth5year_6064 <- NA
d$birth5year_6064[d$birth5year != "1960-64" & !is.na(d$birth5year)] <- 0
d$birth5year_6064[d$birth5year == "1960-64" & !is.na(d$birth5year)] <- 1

d$birth5year_5559 <- NA
d$birth5year_5559[d$birth5year != "1955-59" & !is.na(d$birth5year)] <- 0
d$birth5year_5559[d$birth5year == "1955-59" & !is.na(d$birth5year)] <- 1

d$birth5year_5054 <- NA
d$birth5year_5054[d$birth5year != "1950-54" & !is.na(d$birth5year)] <- 0
d$birth5year_5054[d$birth5year == "1950-54" & !is.na(d$birth5year)] <- 1

d$birth5year_4549 <- NA
d$birth5year_4549[d$birth5year != "1945-49" & !is.na(d$birth5year)] <- 0
d$birth5year_4549[d$birth5year == "1945-49" & !is.na(d$birth5year)] <- 1

d$birth5year_4044 <- NA
d$birth5year_4044[d$birth5year != "1940-44" & !is.na(d$birth5year)] <- 0
d$birth5year_4044[d$birth5year == "1940-44" & !is.na(d$birth5year)] <- 1


### output container
out <- data.frame(matrix(nrow=19, ncol=3))
colnames(out) <- c("Variable", "Mean", "SD")

out$Variable <- c("Democracy satisfaction", "Country's current economic situation", "Own current job situation",
                  "Country's future economic situation", "Own future job situation", 
                  "Born before 1975", "CEE", "Western/Northern Europe",
                  "Born: 1990+", "Born: 1985-89", "Born: 1980-84", "Born: 1975-79", "Born: 1970-74",
                  "Born: 1965-69", "Born: 1960-64", "Born: 1955-59", "Born: 1950-54", 
                  "Born: 1945-49", "Born: 1940-44")


### implement function / get descriptive statistics
out[1,2:3] <- desc(dv = "demsatis")

out[2,2:3] <- desc(dv = "econsit_country")
out[3,2:3] <- desc(dv = "jobsit_personal")
out[4,2:3] <- desc(dv = "expect_econsit")
out[5,2:3] <- desc(dv = "expect_jobsit")

out[6,2:3] <- desc(dv = "born_before_1975")

out[7,2:3] <- desc(dv = "region_CEE")
out[8,2:3] <- desc(dv = "region_West")

out[9,2:3] <- desc(dv = "birth5year_1990")
out[10,2:3] <- desc(dv = "birth5year_8589")
out[11,2:3] <- desc(dv = "birth5year_8084")
out[12,2:3] <- desc(dv = "birth5year_7579")
out[13,2:3] <- desc(dv = "birth5year_7074")
out[14,2:3] <- desc(dv = "birth5year_6569")
out[15,2:3] <- desc(dv = "birth5year_6064")
out[16,2:3] <- desc(dv = "birth5year_5559")
out[17,2:3] <- desc(dv = "birth5year_5054")
out[18,2:3] <- desc(dv = "birth5year_4549")
out[19,2:3] <- desc(dv = "birth5year_4044")

### output
print(xtable(out, digits=3, caption="Descriptive statistics"), 
      caption.placement = "top", include.rownames=FALSE)





# Figure A3 ---------------------------------------------------------------

### re-arrange countries such that they are sorted by low/high UER
uer$country2 <- NA
uer$country2[uer$country == "Czechia"]  <- 1
uer$country2[uer$country == "Estonia"]  <- 2
uer$country2[uer$country == "Hungary"]  <- 3
uer$country2[uer$country == "Romania"]  <- 4
uer$country2[uer$country == "Slovenia"]  <- 5

uer$country2[uer$country == "Bulgaria"]  <- 6
uer$country2[uer$country == "Latvia"]  <- 7
uer$country2[uer$country == "Lithuania"]  <- 8
uer$country2[uer$country == "Poland"]  <- 9
uer$country2[uer$country == "Slovak Republic"]  <- 10

uer$country2 <- factor(uer$country2, levels = 1:10,
                       labels = c("Czechia", "Estonia", "Hungary", "Romania", "Slovenia",
                                  "Bulgaria", "Latvia", "Lithuania", "Poland", "Slovakia"))

### Plot
ggplot() + 
  geom_line(data=uer, 
            aes(x=year, y=unemprate),
            size=1) + 
  facet_wrap(~country2, ncol=5) +
  scale_x_continuous(limits=c(1990,2005), breaks=seq(1990,2002,4)) +
  scale_y_continuous(limits = c(0,22)) +
  xlab("") +
  ylab("Unemployment rate\n(% of labor force)") + 
  theme_classic() +
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=11),
        strip.text = element_text(size=14, face="bold"))





# Figure A4 ---------------------------------------------------------------

### Regressions
# satisfaction with country's current economic situation
m1 <- felm(as.numeric(demsatis) ~ as.numeric(econsit_country)*as.factor(region)*as.factor(birth5year) + 
             as.numeric(econsit_country)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# satisfaction with own current job situation
m2 <- felm(as.numeric(demsatis) ~ as.numeric(jobsit_personal)*as.factor(region)*as.factor(birth5year) + 
             as.numeric(jobsit_personal)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# expectations about country's future economic situation
m3 <- felm(as.numeric(demsatis) ~ as.numeric(expect_econsit)*as.factor(region)*as.factor(birth5year) + 
             as.numeric(expect_econsit)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)

# expectations about own future job situation
m4 <- felm(as.numeric(demsatis) ~ as.numeric(expect_jobsit)*as.factor(region)*as.factor(birth5year) + 
             as.numeric(expect_jobsit)*as.factor(region)*as.factor(agegroup) + as.factor(female) + as.factor(employment) +
             as.factor(marstat) + as.factor(age_education) | country + year | 0 | ID, data=d)


### get coefficients and SEs
# prepare data frame
out <- data.frame(matrix(nrow=44, ncol=4))
colnames(out) <- c("variable", "years", "b", "se")
out$variable <- rep(1:4, each=11)
out$variable <- factor(out$variable, levels=1:4,
                       labels = c("Country's current economic situation",
                                  "Own current job situation",
                                  "Country's future economic situation",
                                  "Own future job situation"))
out$years <- rep(1:11, 4)
out$years <- factor(out$years, levels=1:11,
                    labels =c("1990+", "1985-89", "1980-84", "1975-79", 
                              "1970-74", "1965-69", "1960-64", "1955-59", 
                              "1950-54", "1945-49", "1940-44"))
out$b[out$years == "1990+"] <- 0
out$se[out$years == "1990+"] <- 0

# extract coefs and ses
out[2,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1985-89",1:2]
out[3,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1980-84",1:2]
out[4,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1975-79",1:2]
out[5,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1970-74",1:2]
out[6,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1965-69",1:2]
out[7,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1960-64",1:2]
out[8,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1955-59",1:2]
out[9,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1950-54",1:2]
out[10,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1945-49",1:2]
out[11,3:4] <- summary(m1)$coef["as.numeric(econsit_country):as.factor(region)Eastern Europe:as.factor(birth5year)1940-44",1:2]

out[13,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1985-89",1:2]
out[14,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1980-84",1:2]
out[15,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1975-79",1:2]
out[16,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1970-74",1:2]
out[17,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1965-69",1:2]
out[18,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1960-64",1:2]
out[19,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1955-59",1:2]
out[20,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1950-54",1:2]
out[21,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1945-49",1:2]
out[22,3:4] <- summary(m2)$coef["as.numeric(jobsit_personal):as.factor(region)Eastern Europe:as.factor(birth5year)1940-44",1:2]

out[24,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1985-89",1:2]
out[25,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1980-84",1:2]
out[26,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1975-79",1:2]
out[27,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1970-74",1:2]
out[28,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1965-69",1:2]
out[29,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1960-64",1:2]
out[30,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1955-59",1:2]
out[31,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1950-54",1:2]
out[32,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1945-49",1:2]
out[33,3:4] <- summary(m3)$coef["as.numeric(expect_econsit):as.factor(region)Eastern Europe:as.factor(birth5year)1940-44",1:2]

out[35,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1985-89",1:2]
out[36,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1980-84",1:2]
out[37,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1975-79",1:2]
out[38,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1970-74",1:2]
out[39,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1965-69",1:2]
out[40,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1960-64",1:2]
out[41,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1955-59",1:2]
out[42,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1950-54",1:2]
out[43,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1945-49",1:2]
out[44,3:4] <- summary(m4)$coef["as.numeric(expect_jobsit):as.factor(region)Eastern Europe:as.factor(birth5year)1940-44",1:2]



### Plot
ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=years, y=b, shape=variable), 
             position=position_dodge(width=1), size=3) + 
  geom_errorbar(data=out, aes(x=years, ymin=b-1.96*se, ymax=b+1.96*se),
                width=0, position=position_dodge(width=1), size=1) +
  xlab("Birth cohort") +
  ylab("Estimated difference in the association\nbetween economic evaluations and democracy satisfaction\nbetween Eastern and Western Europeans\nby birth cohort") +
  facet_wrap(~variable, scales="free_y") +
  theme_classic() + 
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=12, angle=90, hjust=1, vjust=.5),
        axis.text.y = element_text(size=12),
        strip.text = element_text(size=14, face="bold"),
        legend.position = "none")




